

% Define density fzeta
fzeta = @(zzeta) ...
            ( zzeta <= Params.zetaU ) .* ( zzeta >= Params.zetaL ) ...
         .* ( ( Params.TailShiftZeta + ...
                ( ( zzeta - Params.zetaL ) / ( Params.zetaU - Params.zetaL ) ).^( Params.p1zeta - 1 ) ...
                .* ( ( Params.zetaU - zzeta ) / ( Params.zetaU - Params.zetaL ) ).^( Params.p2zeta - 1 ) ...
                ) ...
                ./ ( Params.TailShiftZeta + gamma( Params.p1zeta ) * gamma( Params.p2zeta ) / gamma( Params.p1zeta + Params.p2zeta ) ) ...
             ) ...
            + ( ( zzeta > Params.zetaU ) | ( zzeta < Params.zetaL ) ) * Params.TailShiftZeta ;
